! (C) Copyright 2005- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

PROGRAM CREATE_GRID_INFO
USE NETCDF
USE GRIB_API
IMPLICIT NONE

CHARACTER(LEN=500) :: LSM_FILE,NETCDF_FILE,CLAKE_FILE
INTEGER            :: IGRIB,IFILE,NCID,VARID,NPOINTS,TNPOINTS,NLAT,JP,JPP,NLAKE
REAL*8,DIMENSION(:),ALLOCATABLE :: TLAT,TLON,LAT,LON,VALUES,LATS,VALUES2
REAL*8,DIMENSION(:,:),ALLOCATABLE :: XCORNER,YCORNER
INTEGER,DIMENSION(:),ALLOCATABLE :: NLONS,POINTS,LPOINTS
REAL*8             :: LSM_MIN,LAT1,LATL,DELLAT,DELLON

NAMELIST /INPUT/LSM_FILE,CLAKE_FILE,NETCDF_FILE

OPEN(10,FILE='input.nam')
READ(10,INPUT)
CLOSE(10)

!LSM_FILE='lsm_N128.grib'
!CLAKE_FILE='clake.grib'
!NETCDF_FILE='infoN128.nc'
LSM_MIN=0.5


IFILE=2
PRINT*,'OPENING LSM FILE ',TRIM(LSM_FILE)
CALL GRIB_OPEN_FILE(IFILE,TRIM(LSM_FILE),'R')
CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB)

CALL GRIB_GET(IGRIB,'numberOfPoints',TNPOINTS)
ALLOCATE(TLAT(TNPOINTS))
ALLOCATE(TLON(TNPOINTS))
ALLOCATE(VALUES(TNPOINTS))

CALL GRIB_GET_DATA(IGRIB,TLAT,TLON,VALUES)
CALL GRIB_GET(IGRIB,'numberOfPointsAlongAMeridian',NLAT)
CALL GRIB_GET(IGRIB,'latitudeOfFirstGridPointInDegrees',LAT1)
CALL GRIB_GET(IGRIB,'latitudeOfLastGridPointInDegrees',LATL)
DELLAT=(LAT1-LATL)/(NLAT-1.)
ALLOCATE(LATS(NLAT))
ALLOCATE(NLONS(NLAT))
CALL GRIB_GET(IGRIB,'pl',NLONS)
NPOINTS=COUNT(VALUES>LSM_MIN)
ALLOCATE(POINTS(NPOINTS))
ALLOCATE(LAT(NPOINTS))
ALLOCATE(LON(NPOINTS))
ALLOCATE(XCORNER(TNPOINTS,4))
ALLOCATE(YCORNER(TNPOINTS,4))
JPP=0
DO JP=1,TNPOINTS
  IF( VALUES(JP) > LSM_MIN ) THEN
    JPP=JPP+1
    POINTS(JPP)=JP
    LAT(JPP)=TLAT(JP)
    LON(JPP)=TLON(JP)
  ENDIF
ENDDO
DO JP=1,NLAT
  LATS(JP)=LAT1-(JP-1)*DELLAT
!   print*,JP,LATS(JP),NLONS(JP)
ENDDO

JPP=1
DO JP=1,TNPOINTS
  IF ( JP > SUM(NLONS(1:JPP)) ) JPP=JPP+1
  YCORNER(JP,1)=TLAT(JP)-DELLAT*0.5
  YCORNER(JP,2)=TLAT(JP)-DELLAT*0.5
  YCORNER(JP,3)=TLAT(JP)+DELLAT*0.5
  YCORNER(JP,4)=TLAT(JP)+DELLAT*0.5
  DELLON=360./NLONS(JPP)
  XCORNER(JP,1)=TLON(JP)-DELLON*0.5
  XCORNER(JP,2)=TLON(JP)+DELLON*0.5
  XCORNER(JP,3)=TLON(JP)+DELLON*0.5
  XCORNER(JP,4)=TLON(JP)-DELLON*0.5
!   PRINT*,JP,TLAT(JP),TLON(JP),NLONS(JPP)
ENDDO

CALL GRIB_CLOSE_FILE(IFILE)

IFILE=3
PRINT*,'OPENING CLAKE LSM FILE ',TRIM(CLAKE_FILE)
  CALL GRIB_OPEN_FILE(IFILE,TRIM(CLAKE_FILE),'R')
CALL GRIB_NEW_FROM_FILE(IFILE,IGRIB)
ALLOCATE(VALUES2(TNPOINTS))
CALL GRIB_GET(IGRIB,'values',VALUES2)
NLAKE=COUNT(VALUES2>LSM_MIN)
ALLOCATE(LPOINTS(NLAKE))
CALL GRIB_CLOSE_FILE(IFILE)
JPP=0
DO JP=1,TNPOINTS
  IF ( VALUES2(JP) > LSM_MIN) THEN
    JPP=JPP+1
    LPOINTS(JPP)=JP
   ENDIF
ENDDO



PRINT*,'TOTAL GRID POINTS=',TNPOINTS
PRINT*,'LAND POINTS      =',NPOINTS
PRINT*,'LAKE POINTS      =',NLAKE
PRINT*,'NUMBER LATITUDES =',NLAT
PRINT*,LAT1,LATL,DELLAT

CALL NCERROR( NF90_OPEN(TRIM(NETCDF_FILE),NF90_WRITE,NCID))

CALL NCERROR( NF90_INQ_VARID(NCID,'lon',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,TLON,(/1/),(/TNPOINTS/)) )
CALL NCERROR( NF90_INQ_VARID(NCID,'lat',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,TLAT,(/1/),(/TNPOINTS/)) )
CALL NCERROR( NF90_INQ_VARID(NCID,'lsm',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,VALUES,(/1/),(/TNPOINTS/)) )
CALL NCERROR( NF90_INQ_VARID(NCID,'cl',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,VALUES2,(/1/),(/TNPOINTS/)) )


CALL NCERROR( NF90_INQ_VARID(NCID,'lats',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,LATS,(/1/),(/NLAT/)) )
CALL NCERROR( NF90_INQ_VARID(NCID,'nlons',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,NLONS,(/1/),(/NLAT/)) )

CALL NCERROR( NF90_INQ_VARID(NCID,'xcorners',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,XCORNER,(/1,1/),(/TNPOINTS,4/)) )
CALL NCERROR( NF90_INQ_VARID(NCID,'ycorners',VARID) )
CALL NCERROR( NF90_PUT_VAR(NCID,VARID,YCORNER,(/1,1/),(/TNPOINTS,4/)) )


CALL NCERROR( NF90_CLOSE(NCID) )


CONTAINS 
SUBROUTINE NCERROR(STATUS,STRING)
USE NETCDF
IMPLICIT NONE

INTEGER,INTENT(IN) :: STATUS
CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: STRING 

IF ( STATUS /= 0 ) THEN 
  PRINT *, TRIM(NF90_STRERROR(STATUS))
  IF( PRESENT(STRING) ) PRINT*,TRIM(STRING)
  STOP
ENDIF

END SUBROUTINE NCERROR

END PROGRAM CREATE_GRID_INFO
